rm(list=ls(all=TRUE))
setwd() # input working directory
#install.packages("pacman")
pacman::p_load(dplyr, tidyr, tidyverse, infer, broom,
               descr, stargazer, gtools, boot, InteractionPoweR, pwr)

###### Load & clean data #####

data <- read.csv("Casler JCR data final.csv")

##### Covariates #####

data$peace_strength <- ifelse(data$Q5_1=="Strongly disagree",1,
                       ifelse(data$Q5_1=="Somewhat disagree",2,
                       ifelse(data$Q5_1=="Neither agree nor disagree",3,
                       ifelse(data$Q5_1=="Somewhat agree",4,
                       ifelse(data$Q5_1=="Strongly agree",5,NA)))))

data$strike_first <- ifelse(data$Q5_2=="Strongly disagree",1,
                     ifelse(data$Q5_2=="Somewhat disagree",2,
                     ifelse(data$Q5_2=="Neither agree nor disagree",3,
                     ifelse(data$Q5_2=="Somewhat agree",4,
                     ifelse(data$Q5_2=="Strongly agree",5,NA)))))

data$us_influence <- ifelse(data$Q5_3=="Strongly disagree",1,
                     ifelse(data$Q5_3=="Somewhat disagree",2,
                     ifelse(data$Q5_3=="Neither agree nor disagree",3,
                     ifelse(data$Q5_3=="Somewhat agree",4,
                     ifelse(data$Q5_3=="Strongly agree",5,NA)))))

data$force <- ifelse(data$Q5_4=="Strongly disagree",5,
                     ifelse(data$Q5_4=="Somewhat disagree",4,
                     ifelse(data$Q5_4=="Neither agree nor disagree",3,
                     ifelse(data$Q5_4=="Somewhat agree",2,
                     ifelse(data$Q5_4=="Strongly agree",1,NA)))))

data$assert <- (data$peace_strength + data$strike_first + data$us_influence + data$force)/4

data$un <- ifelse(data$Q42_1=="Strongly disagree",1,
           ifelse(data$Q42_1=="Somewhat disagree",2,
           ifelse(data$Q42_1=="Neither agree nor disagree",3,
           ifelse(data$Q42_1=="Somewhat agree",4,
           ifelse(data$Q42_1=="Strongly agree",5,NA)))))

data$others <- ifelse(data$Q42_2=="Strongly disagree",1,
               ifelse(data$Q42_2=="Somewhat disagree",2,
               ifelse(data$Q42_2=="Neither agree nor disagree",3,
               ifelse(data$Q42_2=="Somewhat agree",4,
               ifelse(data$Q42_2=="Strongly agree",5,NA)))))

data$active <- ifelse(data$Q42_3=="Strongly disagree",1,
               ifelse(data$Q42_3=="Somewhat disagree",2,
               ifelse(data$Q42_3=="Neither agree nor disagree",3,
               ifelse(data$Q42_3=="Somewhat agree",4,
               ifelse(data$Q42_3=="Strongly agree",5,NA)))))

data$home <- ifelse(data$Q5_4=="Strongly disagree",5,
             ifelse(data$Q5_4=="Somewhat disagree",4,
             ifelse(data$Q5_4=="Neither agree nor disagree",3,
             ifelse(data$Q5_4=="Somewhat agree",2,
             ifelse(data$Q5_4=="Strongly agree",1,NA)))))

data$coop <- (data$un + data$others + data$active + data$home)/4

data$current_civFA <- ifelse(data$Q6=="Civilian foreign affairs agency",1,0)
data$current_civDom <- ifelse(data$Q6=="Civilian domestic affairs agency",1,0)
data$current_civDef <- ifelse(data$Q6=="Civilian defense agency",1,0)
data$current_mil <- ifelse(data$Q6=="Armed services",1,0)
data$current_pol <- ifelse(data$Q6=="Political party",1,0)
data$current_priv <- ifelse(data$Q6=="Private sector",1,0)
data$current_other <- ifelse(data$Q6=="Other",1,0)

data$prev_civFA <- ifelse(grepl("Civilian foreign affairs agency", data$Q7, fixed = T),1,0)
data$prev_civDom <- ifelse(grepl("Civilian domestic affairs agency", data$Q7, fixed = T),1,0)
data$prev_civDef <- ifelse(grepl("Civilian defense agency", data$Q7, fixed = T),1,0)
data$prev_mil <- ifelse(grepl("Armed services", data$Q7, fixed = T),1,0)
data$prev_pol <- ifelse(grepl("Political party", data$Q7, fixed = T),1,0)
data$prev_priv <- ifelse(grepl("Private sector", data$Q7, fixed = T),1,0)
data$prev_other <- ifelse(grepl("Other", data$Q7, fixed = T),1,0)

data$civFA <- ifelse(data$current_civFA==1 | data$prev_civFA==1,1,0)
data$civDom <- ifelse(data$current_civDom==1 | data$prev_civDom==1,1,0)
data$civDef <- ifelse(data$current_civDef==1 | data$prev_civDef==1,1,0)
data$mil <- ifelse(data$current_mil==1 | data$prev_mil==1,1,0)
data$pol <- ifelse(data$current_pol==1 | data$prev_pol==1,1,0)
data$priv <- ifelse(data$current_priv==1 | data$prev_priv==1,1,0)
data$other <- ifelse(data$current_other==1 | data$prev_other==1,1,0)

data$MilExp <- ifelse(data$Q28=="",0,1)

data$dod <- ifelse(data$MilExp==1 | data$current_civDef==1 | data$prev_civDef==1,1,0)

data$MilYears <- data$Q8
data$CivYears <- data$Q17

data$emp.current_ops <- ifelse(grepl("Operations", data$Q13, fixed = T),1,0)
data$emp.current_strat <- ifelse(grepl("Strategy", data$Q13, fixed = T),1,0)
data$emp.current_policy <- ifelse(grepl("Policy", data$Q13, fixed = T),1,0)
data$emp.current_admin <- ifelse(grepl("Administration", data$Q13, fixed = T),1,0)
data$emp.current_other <- ifelse(grepl("Operations", data$Q13, fixed = T),1,0)

data$emp.prev_ops <- ifelse(grepl("Operations", data$Q15, fixed = T),1,0)
data$emp.prev_strat <- ifelse(grepl("Strategy", data$Q15, fixed = T),1,0)
data$emp.prev_policy <- ifelse(grepl("Policy", data$Q15, fixed = T),1,0)
data$emp.prev_admin <- ifelse(grepl("Administration", data$Q15, fixed = T),1,0)
data$emp.prev_other <- ifelse(grepl("Operations", data$Q15, fixed = T),1,0)

data$ops <- ifelse(data$emp.current_ops==1 | data$emp.prev_ops==1,1,0)
data$strat <- ifelse(data$emp.current_strat==1 | data$emp.prev_strat==1,1,0)
data$policy <- ifelse(data$emp.current_policy==1 | data$emp.prev_policy==1,1,0)
data$admin <- ifelse(data$emp.current_admin==1 | data$emp.prev_admin==1,1,0)
data$other <- ifelse(data$emp.current_ops==1 | data$emp.prev_ops==1,1,0)

data$enlisted <- ifelse(data$Q28=="E-1",1,
                 ifelse(data$Q28=="E-2",2,
                 ifelse(data$Q28=="E-3",3,
                 ifelse(data$Q28=="E-4",4,
                 ifelse(data$Q28=="E-5",5,
                 ifelse(data$Q28=="E-6",6,
                 ifelse(data$Q28=="E-7",7,
                 ifelse(data$Q28=="E-8",8,
                 ifelse(data$Q28=="E-9",9,0)))))))))
data$enlisted_bin <- ifelse(data$enlisted!=0,1,0)

data$warrant <- ifelse(data$Q28=="W-1",1,
                 ifelse(data$Q28=="W-2",2,
                 ifelse(data$Q28=="W-3",3,
                 ifelse(data$Q28=="W-4",4,
                 ifelse(data$Q28=="W-5",5,0)))))
data$warrant_bin <- ifelse(data$warrant!=0,1,0)

data$officer <- ifelse(data$Q28=="O-1",1,
                ifelse(data$Q28=="O-2",2,
                ifelse(data$Q28=="O-3",3,
                ifelse(data$Q28=="O-4",4,
                ifelse(data$Q28=="O-5",5,
                ifelse(data$Q28=="O-6",6,
                ifelse(data$Q28=="O-7",7,
                ifelse(data$Q28=="O-8",8,
                ifelse(data$Q28=="O-9",9,0)))))))))
data$officer_bin <- ifelse(data$officer!=0,1,0)

data$Combat <- ifelse(data$Q41=="Yes",1,0)

data$gender <- ifelse(data$Q2=="Female",1,0)

data$educ <- ifelse(data$Q3=="No formal education",1,
             ifelse(data$Q3=="Less than a complete high school education",2,
             ifelse(data$Q3=="Complete high school education",3,
             ifelse(data$Q3=="Some university-level education, without degree",4,
             ifelse(data$Q3=="University-level education, with degree",5,
             ifelse(data$Q3=="Some post-graduate education, without degree",6,
             ifelse(data$Q3=="Post-graduate education, with degree",7,NA)))))))

data$age <- as.numeric(data$Q40)

data$citizen <- ifelse(data$Q4=="Yes",1,0)

###### Treatments #####

data$Reputation <- ifelse(data$FL_17_DO=="Reputation",1,0) # reputation to uphold but no capabilities
data$ReputationPlusCapabilities <- ifelse(data$FL_17_DO=="ReputationPlusCapabilities",1,0) # reputation to uphold and capabilities
data$Control <- ifelse(data$FL_17_DO=="Control",1,0) # no reputation to uphold and no capabilities -- baseline
data$Capabilities <- ifelse(data$FL_17_DO=="Capabilities",1,0) # no reputation to uphold and capabilities 

data$RepCombined <- ifelse(data$Reputation==1 | data$ReputationPlusCapabilities==1,1,0)
data$RepControl <- ifelse(data$Capabilities==1 | data$Control==1,1,0)
data$CapCombined <- ifelse(data$Capabilities==1 | data$ReputationPlusCapabilities==1,1,0)
data$CapControl <- ifelse(data$Reputation==1 | data$Control==1,1,0)

data$treatC <- ifelse(data$CapControl==1,0,
               ifelse(data$CapCombined==1,1,NA))

data$treatR <- ifelse(data$RepControl==1,0,
               ifelse(data$RepCombined==1,1,NA))

###### Outcomes #####

data$dv_force <- ifelse(data$Q41.1=="Strongly oppose",1,
                 ifelse(data$Q41.1=="Somewhat oppose",2,
                 ifelse(data$Q41.1=="Neither support nor oppose",3,
                 ifelse(data$Q41.1=="Somewhat support",4,
                 ifelse(data$Q41.1=="Strongly support",5,NA)))))

data$dv_marines <- ifelse(data$Q43_1=="Strongly oppose",1,
                   ifelse(data$Q43_1=="Somewhat oppose",2,
                   ifelse(data$Q43_1=="Neither support nor oppose",3,
                   ifelse(data$Q43_1=="Somewhat support",4,
                   ifelse(data$Q43_1=="Strongly support",5,NA)))))

data$dv_flyover <- ifelse(data$Q43_2=="Strongly oppose",1,
                   ifelse(data$Q43_2=="Somewhat oppose",2,
                   ifelse(data$Q43_2=="Neither support nor oppose",3,
                   ifelse(data$Q43_2=="Somewhat support",4,
                   ifelse(data$Q43_2=="Strongly support",5,NA)))))

data$dv_sanctions <- ifelse(data$Q43_3=="Strongly oppose",1,
                     ifelse(data$Q43_3=="Somewhat oppose",2,
                     ifelse(data$Q43_3=="Neither support nor oppose",3,
                     ifelse(data$Q43_3=="Somewhat support",4,
                     ifelse(data$Q43_3=="Strongly support",5,NA)))))

data$dv_intel <- ifelse(data$Q43_4=="Strongly oppose",1,
                 ifelse(data$Q43_4=="Somewhat oppose",2,
                 ifelse(data$Q43_4=="Neither support nor oppose",3,
                 ifelse(data$Q43_4=="Somewhat support",4,
                 ifelse(data$Q43_4=="Strongly support",5,NA)))))

data$dv_nothing <- ifelse(data$Q43_5=="Strongly oppose",1,
                   ifelse(data$Q43_5=="Somewhat oppose",2,
                   ifelse(data$Q43_5=="Neither support nor oppose",3,
                   ifelse(data$Q43_5=="Somewhat support",4,
                   ifelse(data$Q43_5=="Strongly support",5,NA)))))

###### Moderators #####

data$USrep_med <- ifelse(data$Q58_3=="Strongly disagree",1,
                  ifelse(data$Q58_3=="Somewhat disagree",2,
                  ifelse(data$Q58_3=="Neither agree nor disagree",3,
                  ifelse(data$Q58_3=="Somewhat agree",4,
                  ifelse(data$Q58_3=="Strongly agree",5,NA)))))

data$rep_region <- ifelse(data$Q59_1=="Strongly disagree",1,
                   ifelse(data$Q59_1=="Somewhat disagree",2,
                   ifelse(data$Q59_1=="Neither agree nor disagree",3,
                   ifelse(data$Q59_1=="Somewhat agree",4,
                   ifelse(data$Q59_1=="Strongly agree",5,NA)))))

data$rep_world <- ifelse(data$Q59_2=="Strongly disagree",1,
                  ifelse(data$Q59_2=="Somewhat disagree",2,
                  ifelse(data$Q59_2=="Neither agree nor disagree",3,
                  ifelse(data$Q59_2=="Somewhat agree",4,
                  ifelse(data$Q59_2=="Strongly agree",5,NA)))))

###### Manipulation Checks #####

data$manip <- ifelse(data$Q37=="Westria",1,0)
data$manip2 <- ifelse(data$Q37=="I don't remember" & data$Q61=="Westria",1,0)

#data <- data[data$manip==0,]
#data <- data[data$manip2==0,]

##### Conditional Average Treatment Effects #####

# Write function for calculating bootstrapped quantities of interest
bootTreat2 <- function(dat, dv="dv", B=1500){
  bootResults <- matrix(NA, nrow=B, ncol=6)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    dep.var <- eval(parse(text=paste("temp",dv,sep="$")))
    A <- mean(dep.var[which(temp$RepControl==1)], na.rm=TRUE)
    B <- mean(dep.var[which(temp$RepCombined==1)], na.rm=TRUE)
    C <- mean(dep.var[which(temp$CapControl==1)], na.rm=TRUE)
    D <- mean(dep.var[which(temp$CapCombined==1)], na.rm=TRUE)
    bootResults[i,1] <- (B-A)
    bootResults[i,2] <- (D-C)
    drop(list())
  }
  return(list(dv=dv, boot=bootResults, 
              rep=mean(bootResults[,1], na.rm=TRUE), 
              cap=mean(bootResults[,2], na.rm=TRUE)))
}
set.seed(1234)


##### CATEs Among Diplomats (Table 3) ####

dip <- data[data$MilExp==0 & data$civFA==1,]
dip_desc <- dip[!is.na(dip$dv_force),]

# Force Among Diplomats
mean(na.omit(dip_desc$dv_force[dip_desc$RepControl==1])) # 3.376344
mean(na.omit(dip_desc$dv_force[dip_desc$RepCombined==1])) # 3.684932

mean(na.omit(dip_desc$dv_force[dip_desc$CapControl==1])) # 3.45977
mean(na.omit(dip_desc$dv_force[dip_desc$CapCombined==1])) # 3.56962

set.seed(1234)
prep.full.dip.force <- bootTreat2(dip_desc, dv="dv_force")
prep.full.dip.force$rep # 0.3051043
prep.full.dip.force$cap # 0.1063482

#To calculate p-values:
(1 - length(which(prep.full.dip.force$boot[,1] > 0))/nrow(prep.full.dip.force$boot))*2 # p = 0.02666667
(1 - length(which(prep.full.dip.force$boot[,2] > 0))/nrow(prep.full.dip.force$boot))*2 # p = 0.5013333

# Marines Among Diplomats
mean(na.omit(dip_desc$dv_marines[dip_desc$RepControl==1])) # 3.010753
mean(na.omit(dip_desc$dv_marines[dip_desc$RepCombined==1])) # 3.410959

mean(na.omit(dip_desc$dv_marines[dip_desc$CapControl==1])) # 3.091954
mean(na.omit(dip_desc$dv_marines[dip_desc$CapCombined==1])) # 3.291139

set.seed(1234)
prep.full.dip.marines <- bootTreat2(dip_desc, dv="dv_marines")
prep.full.dip.marines$rep # 0.4013617
prep.full.dip.marines$cap # 0.1968299

#To calculate p-values:
(1 - length(which(prep.full.dip.marines$boot[,1] > 0))/nrow(prep.full.dip.marines$boot))*2 # p = 0.02266667
(1 - length(which(prep.full.dip.marines$boot[,2] > 0))/nrow(prep.full.dip.marines$boot))*2 # p = 0.3093333

# Flyover Among Diplomats
mean(na.omit(dip_desc$dv_flyover[dip_desc$RepControl==1])) # 3.763441
mean(na.omit(dip_desc$dv_flyover[dip_desc$RepCombined==1])) # 3.808219

mean(na.omit(dip_desc$dv_flyover[dip_desc$CapControl==1])) # 3.747126
mean(na.omit(dip_desc$dv_flyover[dip_desc$CapCombined==1])) # 3.822785

set.seed(1234)
prep.full.dip.flyover <- bootTreat2(dip_desc, dv="dv_flyover")
prep.full.dip.flyover$rep # 0.03609195
prep.full.dip.flyover$cap # 0.07488388

#To calculate p-values:
(1 - length(which(prep.full.dip.flyover$boot[,1] > 0))/nrow(prep.full.dip.flyover$boot))*2 # p = 0.8613333
(1 - length(which(prep.full.dip.flyover$boot[,2] > 0))/nrow(prep.full.dip.flyover$boot))*2 # p = 0.7266667

# Sanctions Among Diplomats
mean(na.omit(dip_desc$dv_sanctions[dip_desc$RepControl==1])) # 4.16129
mean(na.omit(dip_desc$dv_sanctions[dip_desc$RepCombined==1])) # 4.589041

mean(na.omit(dip_desc$dv_sanctions[dip_desc$CapControl==1])) # 4.425287
mean(na.omit(dip_desc$dv_sanctions[dip_desc$CapCombined==1])) # 4.265823

set.seed(1234)
prep.full.dip.sanctions <- bootTreat2(dip_desc, dv="dv_sanctions")
prep.full.dip.sanctions$rep # 0.4227026
prep.full.dip.sanctions$cap # -0.1581938

#To calculate p-values:
(1 - length(which(prep.full.dip.sanctions$boot[,1] > 0))/nrow(prep.full.dip.sanctions$boot))*2 # p = 0
(1 - length(which(prep.full.dip.sanctions$boot[,2] < 0))/nrow(prep.full.dip.sanctions$boot))*2 # p = 0.256


##### CATEs Among Military (Table 3) #####

mil <- data[data$officer_bin==1 | data$warrant_bin==1,] # just commissioned officers
mil_desc <- mil[!is.na(mil$dv_force),]
mil_desc <- subset(mil_desc, !(MilYears %in% 0 & Q28 %in% c("O-8","O-6"))) # manually remove two respondents who are career civilians but answered "yes" to armed services experience

# Force Among Military
mean(na.omit(mil_desc$dv_force[mil_desc$RepControl==1])) # 3.876923
mean(na.omit(mil_desc$dv_force[mil_desc$RepCombined==1])) # 3.964912

mean(na.omit(mil_desc$dv_force[mil_desc$CapControl==1])) # 3.87931
mean(na.omit(mil_desc$dv_force[mil_desc$CapCombined==1])) # 3.953125

set.seed(1234)
prep.full.mil.force <- bootTreat2(mil_desc, dv="dv_force")
prep.full.mil.force$rep # 0.08003416
prep.full.mil.force$cap # 0.06867763

#To calculate p-values:
(1 - length(which(prep.full.mil.force$boot[,1] > 0))/nrow(prep.full.mil.force$boot))*2 # p = 0.64
(1 - length(which(prep.full.mil.force$boot[,2] > 0))/nrow(prep.full.mil.force$boot))*2 # p = 0.7173333

# Marines Among Military
mean(na.omit(mil_desc$dv_marines[mil_desc$RepControl==1])) # 3.446154
mean(na.omit(mil_desc$dv_marines[mil_desc$RepCombined==1])) # 3.535714

mean(na.omit(mil_desc$dv_marines[mil_desc$CapControl==1])) # 3.431034
mean(na.omit(mil_desc$dv_marines[mil_desc$CapCombined==1])) # 3.539683

set.seed(1234)
prep.full.mil.marines <- bootTreat2(mil_desc, dv="dv_marines")
prep.full.mil.marines$rep # 0.08556568
prep.full.mil.marines$cap # 0.1015184

#To calculate p-values:
(1 - length(which(prep.full.mil.marines$boot[,1] > 0))/nrow(prep.full.mil.marines$boot))*2 # p = 0.696
(1 - length(which(prep.full.mil.marines$boot[,2] > 0))/nrow(prep.full.mil.marines$boot))*2 # p = 0.6546667

# Flyover Among Military
mean(na.omit(mil_desc$dv_flyover[mil_desc$RepControl==1])) # 3.923077
mean(na.omit(mil_desc$dv_flyover[mil_desc$RepCombined==1])) # 3.821429

mean(na.omit(mil_desc$dv_flyover[mil_desc$CapControl==1])) # 3.706897
mean(na.omit(mil_desc$dv_flyover[mil_desc$CapCombined==1])) # 4.031746

set.seed(1234)
prep.full.mil.flyover <- bootTreat2(mil_desc, dv="dv_flyover")
prep.full.mil.flyover$rep # -0.108891
prep.full.mil.flyover$cap # 0.3167341
#To calculate p-values:
(1 - length(which(prep.full.mil.flyover$boot[,1] > 0))/nrow(prep.full.mil.flyover$boot))*2 # p = 0.596
(1 - length(which(prep.full.mil.flyover$boot[,2] > 0))/nrow(prep.full.mil.flyover$boot))*2 # p = 0.1333333


##### Difference in CATEs #####

data_dipmil <- rbind(mil_desc, dip_desc)
data_dipmil$Military <- ifelse(data_dipmil$officer_bin==1 | data_dipmil$warrant_bin==1,1,0)
data_dipmil$Diplomat <- ifelse(data_dipmil$MilExp==0 & data_dipmil$civFA==1,1,0)

# Force (per Table 3, Diplomat CATE = 0.31; Military CATE = 0.08)
t.test((data_dipmil$dv_force[data_dipmil$RepCombined==1 & data_dipmil$Diplomat==1]-data_dipmil$dv_force[data_dipmil$RepControl==1 & data_dipmil$Diplomat==1]),(data_dipmil$dv_force[data_dipmil$RepCombined==1 & data_dipmil$Diplomat==0]-data_dipmil$dv_force[data_dipmil$RepControl==1 & data_dipmil$Diplomat==0])) # p = 0.27

# Marines (per Table 3, Diplomat CATE = 0.40; Military CATE = 0.09)
t.test((data_dipmil$dv_marines[data_dipmil$RepCombined==1 & data_dipmil$Diplomat==1]-data_dipmil$dv_marines[data_dipmil$RepControl==1 & data_dipmil$Diplomat==1]),(data_dipmil$dv_marines[data_dipmil$RepCombined==1 & data_dipmil$Diplomat==0]-data_dipmil$dv_marines[data_dipmil$RepControl==1 & data_dipmil$Diplomat==0])) # p = 0.17

# Flyover (per Table 3, Diplomat CATE = 0.07; Military CATE = 0.32)
t.test((data_dipmil$dv_flyover[data_dipmil$CapCombined==1 & data_dipmil$Diplomat==1]-data_dipmil$dv_flyover[data_dipmil$CapControl==1 & data_dipmil$Diplomat==1]),(data_dipmil$dv_flyover[data_dipmil$CapCombined==1 & data_dipmil$Diplomat==0]-data_dipmil$dv_flyover[data_dipmil$CapControl==1 & data_dipmil$Diplomat==0])) # p = 0.29


##### Pooled Mean Support for Policy Options (Table 4) #####

# Force
mean(na.omit(data_dipmil$dv_force[data_dipmil$Control==1])) # 3.626506
mean(na.omit(data_dipmil$dv_force[data_dipmil$Capabilities==1])) # 3.533333
mean(na.omit(data_dipmil$dv_force[data_dipmil$Reputation==1])) # 3.629032
mean(na.omit(data_dipmil$dv_force[data_dipmil$ReputationPlusCapabilities==1])) # 3.970588
mean(na.omit(data_dipmil$dv_force[data_dipmil$ReputationPlusCapabilities==1]))-mean(na.omit(data_dipmil$dv_force[data_dipmil$Control==1])) # 0.3440822
t.test(data_dipmil$dv_force[data_dipmil$ReputationPlusCapabilities == 1], data_dipmil$dv_force[data_dipmil$Control == 1]) # p = 0.03303

# Marines
mean(na.omit(data_dipmil$dv_marines[data_dipmil$Control==1])) # 3.228916
mean(na.omit(data_dipmil$dv_marines[data_dipmil$Capabilities==1])) # 3.146667
mean(na.omit(data_dipmil$dv_marines[data_dipmil$Reputation==1])) # 3.225806
mean(na.omit(data_dipmil$dv_marines[data_dipmil$ReputationPlusCapabilities==1])) # 3.686567
mean(na.omit(data_dipmil$dv_marines[data_dipmil$ReputationPlusCapabilities==1]))-mean(na.omit(data_dipmil$dv_marines[data_dipmil$Control==1])) # 0.4576515
t.test(data_dipmil$dv_marines[data_dipmil$ReputationPlusCapabilities == 1], data_dipmil$dv_marines[data_dipmil$Control == 1]) # p = 0.02219

#Flyover
mean(na.omit(data_dipmil$dv_flyover[data_dipmil$Control==1])) # 3.843373
mean(na.omit(data_dipmil$dv_flyover[data_dipmil$Capabilities==1])) # 3.813333
mean(na.omit(data_dipmil$dv_flyover[data_dipmil$Reputation==1])) # 3.580645
mean(na.omit(data_dipmil$dv_flyover[data_dipmil$ReputationPlusCapabilities==1])) # 4.043478
mean(na.omit(data_dipmil$dv_flyover[data_dipmil$ReputationPlusCapabilities==1]))-mean(na.omit(data_dipmil$dv_flyover[data_dipmil$Control==1])) # 0.1864773
t.test(data_dipmil$dv_flyover[data_dipmil$ReputationPlusCapabilities == 1], data_dipmil$dv_flyover[data_dipmil$Control == 1]) # p = 0.3362


##### Post-Treatment Items Concerning Reputational Stakes #####

# Main factor: U.S. reputation for reliability
mean(na.omit(dip_desc$USrep_med[dip_desc$RepControl==1])) # 4.022472
mean(na.omit(dip_desc$USrep_med[dip_desc$RepCombined==1])) # 4.138889
t.test(dip_desc$USrep_med[dip_desc$RepCombined == 1], dip_desc$USrep_med[dip_desc$RepControl == 1]) # p = 0.4212

mean(na.omit(mil_desc$USrep_med[mil_desc$RepControl==1])) # 4.171875
mean(na.omit(mil_desc$USrep_med[mil_desc$RepCombined==1])) # 3.982143
t.test(mil_desc$USrep_med[mil_desc$RepCombined == 1], mil_desc$USrep_med[mil_desc$RepControl == 1]) # p = 0.3554

# Regional reputation
mean(na.omit(dip_desc$rep_region[dip_desc$RepControl==1])) # 4
mean(na.omit(dip_desc$rep_region[dip_desc$RepCombined==1])) # 4.194444
t.test(dip_desc$rep_region[dip_desc$RepCombined == 1], dip_desc$rep_region[dip_desc$RepControl == 1]) # p = 0.188

mean(na.omit(mil_desc$rep_region[mil_desc$RepControl==1])) # 4.28125
mean(na.omit(mil_desc$rep_region[mil_desc$RepCombined==1])) # 4.436364
t.test(mil_desc$rep_region[mil_desc$RepCombined == 1], mil_desc$rep_region[mil_desc$RepControl == 1]) # p = 0.3194

# World reputation

mean(na.omit(dip_desc$rep_world[dip_desc$RepControl==1])) # 3.58427
mean(na.omit(dip_desc$rep_world[dip_desc$RepCombined==1])) # 3.763889
t.test(dip_desc$rep_world[dip_desc$RepCombined == 1], dip_desc$rep_world[dip_desc$RepControl == 1]) # p = 0.2428

mean(na.omit(mil_desc$rep_world[mil_desc$RepControl==1])) # 4.140625
mean(na.omit(mil_desc$rep_world[mil_desc$RepCombined==1])) # 4.163636
t.test(mil_desc$rep_world[mil_desc$RepCombined == 1], mil_desc$rep_world[mil_desc$RepControl == 1]) # p = 0.8852


##### Descriptive Stats for Diplomats (Tables A1, A8, A9) #####

dip_desc2 <- data.frame(dip_desc$age,dip_desc$gender,dip_desc$educ,dip_desc$assert,dip_desc$coop,dip_desc$CivYears,dip_desc$civDom,dip_desc$civDef,dip_desc$priv,dip_desc$other)
covs_dip <- c("Age","Gender","Education","Militant Assertiveness","Cooperative Internationalism","Years of Dip Experience","Domestic Affairs","Defense Affairs","Private Sector","Other")

# Table A1
stargazer(dip_desc2,
          type = "text",
          omit.summary.stat = c("p25", "p75"),
          style = "ajps",
          covariate.labels = covs_dip) 
# diplomat sample skews male, but less so than the military, and is a bit younger but has the same rough level of education. less hawkish and more internationalist, unsurprisingly. more than 60 percent indicate working at a civilian foreign affairs agency.

# Compute demographics for diplomats -- Table A8
mean(na.omit(dip_desc$age[dip_desc$CapControl==1])) # 42.91954
mean(na.omit(dip_desc$age[dip_desc$CapCombined==1])) # 44.32911
mean(na.omit(dip_desc$age[dip_desc$RepControl==1])) # 41.67742
mean(na.omit(dip_desc$age[dip_desc$RepCombined==1])) # 46.0274

mean(na.omit(dip_desc$educ[dip_desc$CapControl==1])) # 6.517241
mean(na.omit(dip_desc$educ[dip_desc$CapCombined==1])) # 6.392405
mean(na.omit(dip_desc$educ[dip_desc$RepControl==1])) # 6.516129
mean(na.omit(dip_desc$educ[dip_desc$RepCombined==1])) # 6.383562

mean(na.omit(dip_desc$gender[dip_desc$CapControl==1])) # 0.4367816
mean(na.omit(dip_desc$gender[dip_desc$CapCombined==1])) # 0.3797468
mean(na.omit(dip_desc$gender[dip_desc$RepControl==1])) # 0.3763441
mean(na.omit(dip_desc$gender[dip_desc$RepCombined==1])) # 0.4520548

# Perform randomization check for diplomats -- Table A9
dip_balC <- lm(treatC ~ age + educ + gender, data = dip_desc)
dip_balR <- lm(treatR ~ age + educ + gender, data = dip_desc)
stargazer(dip_balC, dip_balR,
          type = "text",
          dep.var.labels = c("Capability","Reputation"),
          covariate.labels = c("Age","Education","Gender"),
          omit.stat = c("adj.rsq","f","ser"))

##### Descriptive Stats for Military (Tables A2, A10, A11) #####

mil_desc2 <- data.frame(mil_desc$age,mil_desc$gender,mil_desc$educ,mil_desc$assert,mil_desc$coop,mil_desc$MilYears,mil_desc$Combat,mil_desc$civFA,mil_desc$civDom,mil_desc$civDef,mil_desc$priv,mil_desc$other)
covs_mil <- c("Age","Gender","Education","Militant Assertiveness","Cooperative Internationalism","Years of Mil. Experience","Combat","Foreign Affairs","Domestic Affairs","Defense Affairs","Private Sector","Other")

# Table A2
stargazer(mil_desc2,
          type = "text",
          omit.summary.stat = c("p25", "p75"),
          style = "ajps",
          covariate.labels = covs_mil) 
# military sample unsurprisingly skews male, highly educated, a bit more hawkish and less internationalist than diplomats, has over 18 years of service

# Compute demographics for military -- Table A10
mean(na.omit(mil_desc$age[mil_desc$CapControl==1])) # 49.65517
mean(na.omit(mil_desc$age[mil_desc$CapCombined==1])) # 51.57812
mean(na.omit(mil_desc$age[mil_desc$RepControl==1])) # 50.98462
mean(na.omit(mil_desc$age[mil_desc$RepCombined==1])) # 50.29825

mean(na.omit(mil_desc$educ[mil_desc$CapControl==1])) # 6.724138
mean(na.omit(mil_desc$educ[mil_desc$CapCombined==1])) # 6.71875
mean(na.omit(mil_desc$educ[mil_desc$RepControl==1])) # 6.723077
mean(na.omit(mil_desc$educ[mil_desc$RepCombined==1])) # 6.719298

mean(na.omit(mil_desc$gender[mil_desc$CapControl==1])) # 0.137931
mean(na.omit(mil_desc$gender[mil_desc$CapCombined==1])) # 0.140625
mean(na.omit(mil_desc$gender[mil_desc$RepControl==1])) # 0.1538462
mean(na.omit(mil_desc$gender[mil_desc$RepCombined==1])) # 0.122807

# Perform randomization check for military -- Table A11
mil_balC <- lm(treatC ~ age + educ + gender, data = mil_desc)
mil_balR <- lm(treatR ~ age + educ + gender, data = mil_desc)
stargazer(mil_balC, mil_balR,
          type = "text",
          dep.var.labels = c("Capability","Reputation"),
          covariate.labels = c("Age","Education","Gender"),
          omit.stat = c("adj.rsq","f","ser"))

##### Descriptive Stats for Full Sample (Tables A3-A7) #####

data_dipmil <- rbind(mil_desc, dip_desc)
data_dipmil$Military <- ifelse(data_dipmil$officer_bin==1 | data_dipmil$warrant_bin==1,1,0)
data_dipmil$Diplomat <- ifelse(data_dipmil$MilExp==0 & data_dipmil$civFA==1,1,0)

sum(data_dipmil$snowball) # count number of responses that did not come from LinkedIn for Table A3
crosstab(data_dipmil$treatC, data_dipmil$Military) # Table A4 (top two rows)
crosstab(data_dipmil$treatR, data_dipmil$Military) # Table A4 (bottom two rows)

sum(data_dipmil$Military)
sum(data_dipmil$Diplomat)

# Compute distribution of current and previous work experience for Table A5
table(data_dipmil$current_civFA)
table(data_dipmil$current_civDom)
table(data_dipmil$current_civDef)
table(data_dipmil$current_mil)
table(data_dipmil$current_pol)
table(data_dipmil$current_priv)
table(data_dipmil$current_other)

table(data_dipmil$prev_civFA)
table(data_dipmil$prev_civDom)
table(data_dipmil$prev_civDef)
table(data_dipmil$prev_mil)
table(data_dipmil$prev_pol)
table(data_dipmil$prev_priv)
table(data_dipmil$prev_other)

# Compute distribution of employment experience and main tasks for Table A6
table(data_dipmil$civFA, data_dipmil$admin)
table(data_dipmil$civFA, data_dipmil$ops)
table(data_dipmil$civFA, data_dipmil$policy)
table(data_dipmil$civFA, data_dipmil$strat)
table(data_dipmil$civFA, data_dipmil$other)

table(data_dipmil$civDom, data_dipmil$admin)
table(data_dipmil$civDom, data_dipmil$ops)
table(data_dipmil$civDom, data_dipmil$policy)
table(data_dipmil$civDom, data_dipmil$strat)
table(data_dipmil$civDom, data_dipmil$other)

table(data_dipmil$civDef, data_dipmil$admin)
table(data_dipmil$civDef, data_dipmil$ops)
table(data_dipmil$civDef, data_dipmil$policy)
table(data_dipmil$civDef, data_dipmil$strat)
table(data_dipmil$civDef, data_dipmil$other)

table(data_dipmil$mil, data_dipmil$admin)
table(data_dipmil$mil, data_dipmil$ops)
table(data_dipmil$mil, data_dipmil$policy)
table(data_dipmil$mil, data_dipmil$strat)
table(data_dipmil$mil, data_dipmil$other)

table(data_dipmil$pol, data_dipmil$admin)
table(data_dipmil$pol, data_dipmil$ops)
table(data_dipmil$pol, data_dipmil$policy)
table(data_dipmil$pol, data_dipmil$strat)
table(data_dipmil$pol, data_dipmil$other)

table(data_dipmil$priv, data_dipmil$admin)
table(data_dipmil$priv, data_dipmil$ops)
table(data_dipmil$priv, data_dipmil$policy)
table(data_dipmil$priv, data_dipmil$strat)
table(data_dipmil$priv, data_dipmil$other)
# the sample isn't just comprised of secretaries or people doing chiefly administrative work

# Compute distribution of most senior official personally briefed for Table A7
table(data_dipmil$Q9)
# many of the folks in the sample have briefed important/senior national leaders 


##### Main Results with Covariate Adjustment (Tables A12-A14) #####

# Create dataframe for DoD robustness
mildod <- data[data$officer_bin==1 | data$warrant_bin==1 | data$civDef==1,] # plus dod
mildod_desc <- mildod[!is.na(mildod$dv_force),]

# Force 
m1a <- lm(dv_force ~ RepCombined + CapCombined + assert + coop + gender + educ + age + CivYears + snowball, data = dip_desc)
#summary(m1a)

m1b <- lm(dv_force ~ RepCombined + CapCombined + assert + coop + gender + educ + age + MilYears + Combat + snowball, data = mil_desc)
#summary(m1b)

m1c <- lm(dv_force ~ RepCombined + CapCombined + assert + coop + gender + educ + age + MilYears + Combat + snowball, data = mildod_desc)
#summary(m1c)

# Marines
m2a <- lm(dv_marines ~ RepCombined + CapCombined + assert + coop + gender + educ + age + CivYears + snowball, data = dip_desc)
#summary(Marines_dip_covs)

m2b <- lm(dv_marines ~ RepCombined + CapCombined + assert + coop + gender + educ + age + MilYears + Combat + snowball, data = mil_desc)
#summary(m2b)

m2c <- lm(dv_marines ~ RepCombined + CapCombined + assert + coop + gender + educ + age + MilYears + Combat + snowball, data = mildod_desc)
#summary(m2c)

# Flyover
m3a <- lm(dv_flyover ~ RepCombined + CapCombined + assert + coop + gender + educ + age + CivYears + snowball, data = dip_desc)
#summary(m3a)

m3b <- lm(dv_flyover ~ RepCombined + CapCombined + assert + coop + gender + educ + age + MilYears + Combat + snowball, data = mil_desc)
#summary(m3b)

m3c <- lm(dv_flyover ~ RepCombined + CapCombined + assert + coop + gender + educ + age + MilYears + Combat + snowball, data = mildod_desc)
#summary(m3c)

# Sanctions
m4a <- lm(dv_sanctions ~ RepCombined + CapCombined + assert + coop + gender + educ + age + CivYears + snowball, data = dip_desc)
#summary(m4a)

m4b <- lm(dv_sanctions ~ RepCombined + CapCombined + assert + coop + gender + educ + age + MilYears + Combat + snowball, data = mil_desc)
#summary(m4b)

m4c <- lm(dv_sanctions ~ RepCombined + CapCombined + assert + coop + gender + educ + age + MilYears + Combat + snowball, data = mildod_desc)
#summary(m4c)

# Intel
m5a <- lm(dv_intel ~ RepCombined + CapCombined + assert + coop + gender + educ + age + CivYears + snowball, data = dip_desc)
#summary(m5a)

m5b <- lm(dv_intel ~ RepCombined + CapCombined + assert + coop + gender + educ + age + MilYears + Combat + snowball, data = mil_desc)
#summary(m5b)

m5c <- lm(dv_intel ~ RepCombined + CapCombined + assert + coop + gender + educ + age + MilYears + Combat + snowball, data = mildod_desc)
#summary(m5c)

# Nothing
m6a <- lm(dv_nothing ~ RepCombined + CapCombined + assert + coop + gender + educ + age + CivYears + snowball, data = dip_desc)
#summary(m6a)

m6b <- lm(dv_nothing ~ RepCombined + CapCombined + assert + coop + gender + educ + age + MilYears + Combat + snowball, data = mil_desc)
#summary(m6b)

m6c <- lm(dv_nothing ~ RepCombined + CapCombined + assert + coop + gender + educ + age + MilYears + Combat + snowball, data = mildod_desc)
#summary(m6c)

# Generate Tables A12-A14

# Policy Support Among Diplomats (Table A12)
stargazer(m1a, m2a, m3a, m4a, m5a, m6a,
          type = "text",
          dep.var.labels = c("Force","Marines","Flyover","Sanctions","Intelligence","Nothing"),
          #column.labels = c("Diplomats","Military"),
          covariate.labels = c("High Reputation","High Capability","Militant Assertiveness","Cooperative Internationalism","Gender","Education","Age","Years in Govt.","Sample"),
          #add.lines = list(c("Controls","Yes","Yes")),
          omit.stat = c("adj.rsq","f","ser"))

# Policy Support Among Military Officers (Table A13)
stargazer(m1b, m2b, m3b, m4b, m5b, m6b,
          type = "text",
          dep.var.labels = c("Force","Marines","Flyover","Sanctions","Intelligence","Nothing"),
          #column.labels = c("Diplomats","Military"),
          covariate.labels = c("High Reputation","High Capability","Militant Assertiveness","Cooperative Internationalism","Gender","Education","Age","Years in Mil.","Combat","Sample"),
          #add.lines = list(c("Controls","Yes","Yes")),
          omit.stat = c("adj.rsq","f","ser"))

# Policy Support Among Military Officers + DOD Civilians
stargazer(m1c, m2c, m3c, m4c, m5c, m6c,
          type = "text",
          dep.var.labels = c("Force","Marines","Flyover","Sanctions","Intelligence","Nothing"),
          #column.labels = c("Diplomats","Military"),
          covariate.labels = c("High Reputation","High Capability","Militant Assertiveness","Cooperative Internationalism","Gender","Education","Age","Years in Mil.","Combat","Sample"),
          #add.lines = list(c("Controls","Yes","Yes")),
          omit.stat = c("adj.rsq","f","ser"))


##### Pooled Mean Support for Policy Options, Weighted (Table A15) #####

# Force
Bmean_mil_force <- mean(na.omit(mil_desc$dv_force[mil_desc$Control==1])) # 
Cmean_mil_force <- mean(na.omit(mil_desc$dv_force[mil_desc$Capabilities==1])) # 
Rmean_mil_force <- mean(na.omit(mil_desc$dv_force[mil_desc$Reputation==1])) # 
RPCmean_mil_force <- mean(na.omit(mil_desc$dv_force[mil_desc$ReputationPlusCapabilities==1])) # 

Bmean_dip_force <- mean(na.omit(dip_desc$dv_force[dip_desc$Control==1])) # 
Cmean_dip_force <- mean(na.omit(dip_desc$dv_force[dip_desc$Capabilities==1])) # 
Rmean_dip_force <- mean(na.omit(dip_desc$dv_force[dip_desc$Reputation==1])) # 
RPCmean_dip_force <- mean(na.omit(dip_desc$dv_force[dip_desc$ReputationPlusCapabilities==1])) # 

# Calculate the weighted means

#crosstab(data_dipmil$FL_17_DO, data_dipmil$Military, plot = F)
m_share <- 122/288
d_share <- 166/288

# Force column of Table A15
m_share*Bmean_mil_force + d_share*Bmean_dip_force
m_share*Cmean_mil_force + d_share*Cmean_dip_force
m_share*Rmean_mil_force + d_share*Rmean_dip_force
m_share*RPCmean_mil_force + d_share*RPCmean_dip_force

# Calculate the weighted ATEs

C_mil_force <- mean(na.omit(mil_desc$dv_force[mil_desc$Capabilities==1]))-mean(na.omit(mil_desc$dv_force[mil_desc$Control==1]))
R_mil_force <- mean(na.omit(mil_desc$dv_force[mil_desc$Reputation==1]))-mean(na.omit(mil_desc$dv_force[mil_desc$Control==1]))
RPC_mil_force <- mean(na.omit(mil_desc$dv_force[mil_desc$ReputationPlusCapabilities==1]))-mean(na.omit(mil_desc$dv_force[mil_desc$Control==1]))

C_dip_force <- mean(na.omit(dip_desc$dv_force[dip_desc$Capabilities==1]))-mean(na.omit(dip_desc$dv_force[dip_desc$Control==1]))
R_dip_force <- mean(na.omit(dip_desc$dv_force[dip_desc$Reputation==1]))-mean(na.omit(dip_desc$dv_force[dip_desc$Control==1]))
RPC_dip_force <- mean(na.omit(dip_desc$dv_force[dip_desc$ReputationPlusCapabilities==1]))-mean(na.omit(dip_desc$dv_force[dip_desc$Control==1]))

m_share*C_mil_force + d_share*C_dip_force # -0.12
m_share*R_mil_force + d_share*R_dip_force # -0.07
m_share*RPC_mil_force + d_share*RPC_dip_force # 0.32

# Marines
Bmean_mil_marines <- mean(na.omit(mil_desc$dv_marines[mil_desc$Control==1])) # 
Cmean_mil_marines <- mean(na.omit(mil_desc$dv_marines[mil_desc$Capabilities==1])) # 
Rmean_mil_marines <- mean(na.omit(mil_desc$dv_marines[mil_desc$Reputation==1])) # 
RPCmean_mil_marines <- mean(na.omit(mil_desc$dv_marines[mil_desc$ReputationPlusCapabilities==1])) # 

Bmean_dip_marines <- mean(na.omit(dip_desc$dv_marines[dip_desc$Control==1])) # 
Cmean_dip_marines <- mean(na.omit(dip_desc$dv_marines[dip_desc$Capabilities==1])) # 
Rmean_dip_marines <- mean(na.omit(dip_desc$dv_marines[dip_desc$Reputation==1])) # 
RPCmean_dip_marines <- mean(na.omit(dip_desc$dv_marines[dip_desc$ReputationPlusCapabilities==1])) # 

# Calculate the weighted means

#crosstab(data_dipmil$FL_17_DO, data_dipmil$Military, plot = F)
m_share <- 122/288
d_share <- 166/288

# Marines column of Table A15
m_share*Bmean_mil_marines + d_share*Bmean_dip_marines
m_share*Cmean_mil_marines + d_share*Cmean_dip_marines
m_share*Rmean_mil_marines + d_share*Rmean_dip_marines
m_share*RPCmean_mil_marines + d_share*RPCmean_dip_marines

# Calculate the weighted ATEs

C_mil_marines <- mean(na.omit(mil_desc$dv_marines[mil_desc$Capabilities==1]))-mean(na.omit(mil_desc$dv_marines[mil_desc$Control==1]))
R_mil_marines <- mean(na.omit(mil_desc$dv_marines[mil_desc$Reputation==1]))-mean(na.omit(mil_desc$dv_marines[mil_desc$Control==1]))
RPC_mil_marines <- mean(na.omit(mil_desc$dv_marines[mil_desc$ReputationPlusCapabilities==1]))-mean(na.omit(mil_desc$dv_marines[mil_desc$Control==1]))

C_dip_marines <- mean(na.omit(dip_desc$dv_marines[dip_desc$Capabilities==1]))-mean(na.omit(dip_desc$dv_marines[dip_desc$Control==1]))
R_dip_marines <- mean(na.omit(dip_desc$dv_marines[dip_desc$Reputation==1]))-mean(na.omit(dip_desc$dv_marines[dip_desc$Control==1]))
RPC_dip_marines <- mean(na.omit(dip_desc$dv_marines[dip_desc$ReputationPlusCapabilities==1]))-mean(na.omit(dip_desc$dv_marines[dip_desc$Control==1]))

m_share*C_mil_marines + d_share*C_dip_marines
m_share*R_mil_marines + d_share*R_dip_marines
m_share*RPC_mil_marines + d_share*RPC_dip_marines

# Flyover
Bmean_mil_flyover <- mean(na.omit(mil_desc$dv_flyover[mil_desc$Control==1])) # 
Cmean_mil_flyover <- mean(na.omit(mil_desc$dv_flyover[mil_desc$Capabilities==1])) # 
Rmean_mil_flyover <- mean(na.omit(mil_desc$dv_flyover[mil_desc$Reputation==1])) # 
RPCmean_mil_flyover <- mean(na.omit(mil_desc$dv_flyover[mil_desc$ReputationPlusCapabilities==1])) # 

Bmean_dip_flyover <- mean(na.omit(dip_desc$dv_flyover[dip_desc$Control==1])) # 
Cmean_dip_flyover <- mean(na.omit(dip_desc$dv_flyover[dip_desc$Capabilities==1])) # 
Rmean_dip_flyover <- mean(na.omit(dip_desc$dv_flyover[dip_desc$Reputation==1])) # 
RPCmean_dip_flyover <- mean(na.omit(dip_desc$dv_flyover[dip_desc$ReputationPlusCapabilities==1])) # 

# Calculate the weighted means

#crosstab(data_dipmil$FL_17_DO, data_dipmil$Military, plot = F)
m_share <- 122/288
d_share <- 166/288

# Flyover column of Table A15
m_share*Bmean_mil_flyover + d_share*Bmean_dip_flyover
m_share*Cmean_mil_flyover + d_share*Cmean_dip_flyover
m_share*Rmean_mil_flyover + d_share*Rmean_dip_flyover
m_share*RPCmean_mil_flyover + d_share*RPCmean_dip_flyover

# Calculate the weighted ATEs

C_mil_flyover <- mean(na.omit(mil_desc$dv_flyover[mil_desc$Capabilities==1]))-mean(na.omit(mil_desc$dv_flyover[mil_desc$Control==1]))
R_mil_flyover <- mean(na.omit(mil_desc$dv_flyover[mil_desc$Reputation==1]))-mean(na.omit(mil_desc$dv_flyover[mil_desc$Control==1]))
RPC_mil_flyover <- mean(na.omit(mil_desc$dv_flyover[mil_desc$ReputationPlusCapabilities==1]))-mean(na.omit(mil_desc$dv_flyover[mil_desc$Control==1]))

C_dip_flyover <- mean(na.omit(dip_desc$dv_flyover[dip_desc$Capabilities==1]))-mean(na.omit(dip_desc$dv_flyover[dip_desc$Control==1]))
R_dip_flyover <- mean(na.omit(dip_desc$dv_flyover[dip_desc$Reputation==1]))-mean(na.omit(dip_desc$dv_flyover[dip_desc$Control==1]))
RPC_dip_flyover <- mean(na.omit(dip_desc$dv_flyover[dip_desc$ReputationPlusCapabilities==1]))-mean(na.omit(dip_desc$dv_flyover[dip_desc$Control==1]))

m_share*C_mil_flyover + d_share*C_dip_flyover
m_share*R_mil_flyover + d_share*R_dip_flyover
m_share*RPC_mil_flyover + d_share*RPC_dip_flyover

###### Power Calculations #####

# Diplomats -- Force has 51% power (Hypothesis 1a)
pwr.2p2n.test(h = 0.31,
            n1 = 93,
            n2 = 73,
            sig.level = 0.05,
            alternative = "two.sided")

# with 80% power, can detect an effect as small as 0.44
pwr.2p2n.test(power = 0.80,
              n1 = 93,
              n2 = 73,
              sig.level = 0.05,
              alternative = "two.sided")

# Diplomats -- Marines has 73% power (Hypothesis 1a)
pwr.2p2n.test(h = 0.40,
              n1 = 93,
              n2 = 73,
              sig.level = 0.05,
              alternative = "two.sided")

# with 80% power, can detect an effect as small as 0.44
pwr.2p2n.test(power = 0.80,
              n1 = 93,
              n2 = 73,
              sig.level = 0.05,
              alternative = "two.sided")

# Military -- Flyover has 42% power (Hypothesis 1b)
pwr.2p2n.test(h = 0.32,
              n1 = 64,
              n2 = 58,
              sig.level = 0.05,
              alternative = "two.sided")

# with 80% power, can detect an effect as small as 0.51
pwr.2p2n.test(power = 0.8,
              n1 = 64,
              n2 = 58,
              sig.level = 0.05,
              alternative = "two.sided")

# Baseline versus Reputation Plus Capabilities for Force has 55% power (Hypothesis 2)
pwr.2p2n.test(h = 0.34,
              n1 = 83, # 32+51
              n2 = 68, # 31+37
              sig.level = 0.05,
              alternative = "two.sided")

# Baseline versus Reputation Plus Capabilities for Marines has 80% power (Hypothesis 2)
pwr.2p2n.test(h = 0.46,
              n1 = 83, # 32+51
              n2 = 68, # 31+37
              sig.level = 0.05,
              alternative = "two.sided")

# with 80% power, can detect an effect as small as 0.33
pwr.2p2n.test(power = 0.8,
              n1 = 83,
              n2 = 68,
              sig.level = 0.05,
              alternative = "two.sided")

# Power for interaction using Force as outcome

test_power <- power_interaction_r2(
  alpha = 0.05,             # alpha, for the power analysis
  N = 288,                  # sample size
  r.x1x2.y = 0.23,           # interaction effect to test 
  r.x1.y = cor(data_dipmil$RepCombined, data_dipmil$dv_force),              # correlation between x1 and y
  r.x2.y = cor(data_dipmil$Diplomat,data_dipmil$dv_force),              # correlation between x2 and y
  r.x1.x2 = cor(data_dipmil$Diplomat,data_dipmil$RepCombined),              # correlation between x1 and x2
)

# similar results obtained using simulation
test_power2 <- power_interaction(
  n.iter = 10000,
  alpha = 0.05,             # alpha, for the power analysis
  N = 288,                  # sample size
  r.x1x2.y = 0.23,           # interaction effect to test 
  r.x1.y = cor(data_dipmil$RepCombined, data_dipmil$dv_force),              # correlation between x1 and y 
  r.x2.y = cor(data_dipmil$Diplomat,data_dipmil$dv_force),              # correlation between x2 and y 
  r.x1.x2 = cor(data_dipmil$Diplomat,data_dipmil$RepCombined),              # correlation between x1 and x2
)


###### Follow-Up Study (Table A16) ######

# Load and clean data
data2 <- read.csv("JCR Follow-Up Data.csv")

data2$Reputation <- ifelse(data2$FL_17_DO=="Reputation",1,0)
data2$Baseline <- ifelse(data2$FL_17_DO=="Baseline",1,0)

# Regional reputation ("rr")
data2$rep_region <- ifelse(data2$Q59_1=="Strongly disagree",1,
                    ifelse(data2$Q59_1=="Somewhat disagree",2,
                    ifelse(data2$Q59_1=="Neither agree nor disagree",3,
                    ifelse(data2$Q59_1=="Somewhat agree",4,
                    ifelse(data2$Q59_1=="Strongly agree",5,NA)))))
mean(data2$rep_region[data2$Reputation==1])
mean(data2$rep_region[data2$Baseline==1])
t.test(data2$rep_region[data2$Reputation==1],data2$rep_region[data2$Baseline==1])

# World reputation ("rw")
data2$rep_world <- ifelse(data2$Q59_2=="Strongly disagree",1,
                   ifelse(data2$Q59_2=="Somewhat disagree",2,
                   ifelse(data2$Q59_2=="Neither agree nor disagree",3,
                   ifelse(data2$Q59_2=="Somewhat agree",4,
                   ifelse(data2$Q59_2=="Strongly agree",5,NA)))))
mean(data2$rep_world[data2$Reputation==1])
mean(data2$rep_world[data2$Baseline==1])
t.test(data2$rep_world[data2$Reputation==1],data2$rep_world[data2$Baseline==1])

# Reputation in specific scenario ("rs")
data2$rep_specific <- as.numeric(data2$Q48)
mean(data2$rep_specific[data2$Reputation==1])
mean(data2$rep_specific[data2$Baseline==1])
t.test(data2$rep_specific[data2$Reputation==1],data2$rep_specific[data2$Baseline==1])

# Reputation in general ("rg")
data2$rep_general <- as.numeric(data2$Q49)
mean(data2$rep_general[data2$Reputation==1])
mean(data2$rep_general[data2$Baseline==1])
t.test(data2$rep_general[data2$Reputation==1],data2$rep_general[data2$Baseline==1])

# Calculate bootstrapped effects (Table A16)
bootTreat3 <- function(dat, dv="dv", B=1500){
  bootResults <- matrix(NA, nrow=B, ncol=6)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    dep.var <- eval(parse(text=paste("temp",dv,sep="$")))
    A <- mean(dep.var[which(temp$Baseline==1)], na.rm=TRUE)
    B <- mean(dep.var[which(temp$Reputation==1)], na.rm=TRUE)
    bootResults[i,1] <- (B-A)
    drop(list())
  }
  return(list(dv=dv, boot=bootResults, 
              rep=mean(bootResults[,1], na.rm=TRUE)))
}
set.seed(1234)

prep.full.rr <- bootTreat3(data2, dv="rep_region")
(1 - length(which(prep.full.rr$boot[,1] > 0))/nrow(prep.full.rr$boot))*2

prep.full.rw <- bootTreat3(data2, dv="rep_world")
(1 - length(which(prep.full.rw$boot[,1] > 0))/nrow(prep.full.rw$boot))*2

prep.full.rs <- bootTreat3(data2, dv="rep_specific")
(1 - length(which(prep.full.rs$boot[,1] > 0))/nrow(prep.full.rs$boot))*2

prep.full.rg <- bootTreat3(data2, dv="rep_general")
(1 - length(which(prep.full.rg$boot[,1] > 0))/nrow(prep.full.rg$boot))*2


